#!/usr/bin/perl
# 2012/02/09 2つのsignificant.bedから共通のpeakを数え、横に並べて表示したファイルを作成する
# 2010/08/19 2つのSignificant.bedから共通のPeakの数を数える

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [first bed file] -b [second bed file] -o [output merged bed file]\n";
}

my %opt;
getopts("a:b:o:", \%opt);
my $FILE1 = $opt{a};
my $FILE2 = $opt{b};
#chr1	149450	149900	149690	24880.104	59.71	0	0
#chr1	365000	365390	365180	10971.672	61.641	0	0
#chr1	445850	446210	446030	8530.467	45.645	0	0
#chr1	456350	456770	456560	14966.053	48.099	0	0
#chr1	702920	703340	703130	11036.8	52.434	0	0
#chr1	717950	718370	718160	17702.709	61.471	0	0

my $DUMMY = "\t"x4;

my $FILE_out = $opt{o};

my @First = &readPeaks($FILE1);
my @Second = &readPeaks($FILE2);

sub readPeaks{
	my ($file) = @_;
	my @array;
	my $fh = IO::File->new($file) or die "cannot open $file: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^#/){
			next;
		}
		my (@data) = split /\t/;
		my @temp = @data[0,1,2,3,4];
		push @array, \@temp;
	}
	$fh->close();
	return @array;
}


#---------------------------------------
# 2つのピークのリストを比較する
#---------------------------------------
my $i = 0;
my $j = 0;
my $common = 0;
my $NUM_first = @First;
my $NUM_second = @Second;

my $out = IO::File->new($FILE_out, 'w') or die "canot write $FILE_out: $!";
while($i < $NUM_first or $j < $NUM_second){
	### どちらかが、既にリストが尽きていた場合は、片方だけ出力
	unless($i < $NUM_first){
		$out->print("$DUMMY\t" , join("\t", @{$Second[$j]}), "\n");
		$j++;
		next;
	}
	unless($j < $NUM_second){
		$out->print(join("\t", @{$First[$i]}), "\t$DUMMY\n");
		$i++;
		next;
	}

	my ($chr1, $start1, $end1, @other1) = @{$First[$i]};
	my ($chr2, $start2, $end2, @other2) = @{$Second[$j]};

	if(
		($start2 - 300 < $start1 and $start1 < $end2 + 300) or
		($start1 - 300 < $start2 and $start2 < $end1 + 300)
	){
		$common++;
		$out->print(join("\t", @{$First[$i]}, @{$Second[$j]}) . "\n");
		$i++;
		$j++;
	}else{
		if($chr1 ne $chr2){
			if($chr1 eq 'I' or ($chr1 eq 'II' and $chr2 eq 'III')){
				$out->print(join("\t", @{$First[$i]}), "\t$DUMMY\n");
				$i++;
			}else{
				$out->print("$DUMMY\t" , join("\t", @{$Second[$j]}), "\n");
				$j++;
			}
		}else{
			if($start1 < $start2){
				$out->print(join("\t", @{$First[$i]}), "\t$DUMMY\n");
				$i++;
			}else{
				$out->print("$DUMMY\t" , join("\t", @{$Second[$j]}), "\n");
				$j++;
			}
		}
	}
}
$out->close();

printf "FILE1: %d, Unique: %d, Share: %.2f\n", $NUM_first, $NUM_first - $common, $common / $NUM_first;
printf "FILE2: %d, Unique: %d, Share: %.2f\n", $NUM_second, $NUM_second - $common, $common / $NUM_second;
printf "Common: %d\n", $common;








